function y = fhan(x1, x2, r, h)
d = r*h;
d0 = h*d;
g = x1+h*x2;
absg = abs(g);
a0 = sqrt(d*d+8*r*absg);

if absg > d0
    a = x2 + (a0-d)/2*sign(g);
else
    a = x2 + g/h;
end

absa = abs(a);

if absa >d
    y = -r*sign(a);
else
    y = -r*a/d;
end